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Abstract. In the framework of the second-quantization representation of the master 
equation governing the irreversible epitaxial growth, exact equations describing the 
evolution of the island densities has been obtained. Their decoupling within a mean 
field-type approximation with the unknown correlation functions replaced by capture 
numbers (CNs) has been used to derive a closed set of rate equations. The latter has 
been compared with the exact equations to obtain rigorous definitions of the CNs. The 
CN that describes the nucleation of dimer islands from two mobile monomers has been 
^ | measured in the exact kinetic Monte Carlo simulations with the use of the rigorous 

O ■ definition. Strong disagreement with the literature values calculated within alternative 

techniques has been found, especially at low surface coverage. Plausible causes for the 
discrepancies are suggested. 

Another important result of the rigorous approach is the rate equation with the 
term describing the monomer diffusion. The equation significantly differs from the 

£f) ' widely used equations of this kind known from literature. Arguments in favour of our 

C") , approach are given. 



^ . PACS numbers: 68.55.A-,81.15.Aa 
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1. Introduction 

Lattice gas models (LGMs) with the gas atoms subject to stochastic kinetics are widely 
used in modelling the coherent growth on the surface [HOE]. The efficiency of the 
approach comes from two sources. First, the coherent LGM configuration can be fully 
characterized by the occupation numbers of a discrete lattice. This greatly reduces the 
phase space of the system which in more refined treatments is continuous and thus 
more difficult to deal with. Second, the stochastic kinetics of the LGMs allow for the 
use of the kinetic Monte Carlo (KMC) technique which in the case of discrete variables 
is capable of simulating the temporal evolution of the system on the experimental time 
scale in contrast to more accurate molecular dynamics methods which are restricted to 
very short times because of the computationally demanding continuous phase [2j [3]. 

In practical surface growth studies, however, even the use of the KMC meets with 
difficulties. For example, one of the most frequently studied quantities is the island 
size distribution [2] whose measurement requires gathering statistics on the number of 
islands of all sizes. Because in many cases the islands consist of several hundreds or even 
thousands of atoms [2], reliable statistics can be gathered only in repeated simulations of 
large systems. Therefore, sometimes additional arguments need be invoked to interpret 
the simulated data because of poor statistics [I] . Besides, there exist situations when the 
physics of the problem requires simulation of prohibitively large systems. For example, 
in the cases when the surface diffusivity is very large in comparison with the deposition 
rate, the KMC simulations of systems smaller then some critical value would lead to 
qualitatively incorrect results due to the finite size effects [3 [6] . Accelerated algorithms 
were developed to treat very large systems [TJ [8] but their implementation in concrete 
cases meets with difficulties [9j E]. 

Some of the difficulties with practical use of the KMC technique are a direct 
consequence of its strength, that is, of its ability to give a very detailed description 
of the system under study, such as, e. g., the precise morphologies of the growing 
islands. Such descriptions require considerable resources and efforts to obtain but are 
superfluous in the cases when only a coarse-grained picture of the growth is needed. For 
example, in [10J it was shown that to determine the parameters characterizing the atomic 
diffusivity on the surface in the sub-coalescence growth regime it is quite sufficient to 
gather experimental information only on the island densities with their morphologies 
being irrelevant. 

Simplification of the modelling in such cases has been the subject of extensive 
studies (see, e. g., review articles [HI |2]). It was found that the island densities can be 
calculated with good accuracy in the framework of the mean field-type rate equations 
(REs) whose accurate solution requires only a fraction of the computational effort needed 
for the corresponding KMC study (see [T2J 031 El ELU El HSJ [T6j and references therein). 
But this efficiency comes at the price of the necessity to independently assess many 
parameters or even functions entering the REs. 

The number and the nature of the parameters depend on the specifics of the 
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model under consideration. In the irreversible growth model that will be discussed 
in the present paper the input to the REs is a large (in theory infinite) set of the 
so-called capture numbers (CNs) that describe the rates with which islands of various 
sizes incorporate diffusing atoms (or mobile monomers). The CNs has been the subject 
of extensive investigation since the introduction of the REs in this problem [13]. In 
numerous studies a variety of functional forms and of the factors influencing the CNs 
were researched, such as the power-law dependence on the island size depending on 
the islands morphology [H], on the total island density (see pZT] HE] and references 
therein), the influence of spatial correlations and capture zones, [19] [201 ED [22], of self- 
consistency conditions and total coverage pj)J[23], etc.. Though considerable insights 
has been gained in these studies, fully satisfactory solution has not yet be found so new 
ideas are being put forward in in the search for accurate REs for the epitaxial growth 
[21 ESI [22]. 

One of the difficulties hampering the development of the theory is that in existing 
approaches the quality of the CNs is assessed either indirectly via the quality of the island 
size distributions they produce or with the use of an empirical procedure developed in 
[20] [2] [25] . A major goal of the present study is to propose a way to resolve this difficulty 
by suggesting a method of a direct measurement of the CNs in the KMC simulations. 
Our proposal is based on the observation that the stochastic KMC process can be exactly 
described by the corresponding master equation (ME) [26] [3]. With the use of the ME 
we can derive exact expressions describing the time evolution of the island densities. 
These expressions can be transformed into a closed set of REs by means of decoupling 
of a mean field kind. The CNs that appear in the decoupling can be compared to the 
exact expressions to obtain their precise meaning. With the use of these expressions the 
CNs will be "measured" in the exact KMC simulations. In the course of the derivation 
we will obtain an evolution equation containing the monomer diffusion term also found 
in the self-consistent approaches to the CNs calculation [15] [23] [27] [28]. Because our 
equation differs from the known ones and arguably is more accurate, we consider it as 
our second major achievement in this study. We hope that in the future this equation 
will allow us to improve the self-consistent theory of the CNs JT5] [23] [27] [28] . 

The paper is organized as follows. In the next section we introduce the second 
quantization representation of the stochastic processes that describe the coherent 
epitaxial growth; in section [3] the exact evolution equations are derived and the REs 
obtained via their decoupling are presented in section H] in particular, in section fl~Tl we 
show how the conventional REs known from literature arise in our approach. In section 
Owe will measure the capture rate for the nucleation of dimers in the so-called i=l model 
of irreversible growth [2] . Because the capture rate for the creation of dimers effectively 
defines the total island density [15], its accurate definition may improve the estimates 
of the surface diffusion on the basis of experimental data [10J. In the concluding section 
[6] we briefly discuss the results obtained and possible farther development. 
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2. Master equation in the second quantization representation 

The quantum formalism in application to stochastic LGMs was introduced in [29j E2 ED] 
and subsequently applied by many authors to various stochastic models (for further 
bibliography see J32j [33]) . The formalism does not introduce any new physics into the 
problem but allows one to efficiently deal with huge matrices that appear in the ME 
approach in the kinetic many-body problems using the powerful tools from the second 
quantization or the quantum field theories. 

For concreteness we will illustrate the formalism using the simple model of 
submonolayer growth previously studied in, e. g., [25]. In this model the substrate 
has the square lattice geometry. It is assumed that each site can be occupied by one 
atom at most. Thus, the system consisting of M sites the number of configurations is 
2 M . 

2.1. The Fock representation of the state space 

Usually in LGMs the atomic configuration is characterized by the set of the occupation 
numbers 

to = o,i}, (i) 

where i is the lattice site index. But it is obvious that any quantity acquiring two values 
can be used. In the quantum formalism the state of the system under consideration 
which we denote by the "ket" vector \a) can be characterized by the product over all 
lattice sites 

h = 0, 1 (2) 



and |1>,= 1 , (3) 





where the first vector corresponds to the empty site and the second vector to the occupied 
one. Now, by introducing corresponding conjugate "bra"-vectors via the transposition 
of vectors ([3]) one can see that the set of state vectors (J2J), besides being complete, is 
also orthonormal. 

Transitions between the states of different occupation in ([3]) are performed with the 
use of the creation and annihilation operators of the hard-core bosons [HH [33] 



and Oi = ] (4) 





In the space of vectors (jSJ) of dimensionality 2 M the matrices have the size 2 M x 2 M . 
Therefore, the action of operators (Hj) should be augmented with unit 2x2 matrices 
acting on the sites different from i. For simplicity in derivations below we omit these 
trivial factors. 
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As is easily verified with the use of the explicit representations above, Oj and a+ 
satisfy all the usual properties of the creation and annihilation operators 

0+10)* = 11)^ a,|l> i = 10)^. (5) 

In particular, the operator whose eigenvalues correspond to the occupation numbers ([T]) 
has the conventional form (for simplicity here and below we will use the same notation 
for the operator and for the occupation numbers) 

(6) 

as is easily seen from (jlj). Despite the no-double-occupancy conditions that follow from 
the definitions (J3J) and are reminiscent of the fermions: 

[af, Oj} = 1 and (af) 2 = a 2 = (7) 

(the braces denote the anticommutator), the operators at different sites commute: 

[a+ aj ].^. = 0. (8) 

Finally, by defining the vacuum state of the system as that of the empty lattice 

io>=ni°>« ( 9 ) 

i 

we can represent the state vectors (j2j) of the system as the Fock space 

i«)=( n a * + )i°>- ( io ) 

\Occupied sites / 

2.2. Master equation in the Fock space 

The stochastic LGM of the surface growth that we consider in the present paper belongs 
to the general category of stochastic models whose kinetics satisfy the ME 
HP 

^ = J2 ( WafiPfi ~ W? a P a ) , (11) 

P 

where P a (t) is the probability of the system to be found at time t in state a and W a p 
is the matrix of transition rates between different states of the system. As was shown 
above, in the case of the growth on the surface in the representation of occupation 
vectors ([3]) the matrix has M binary indices and its size is 2 M x 2 M . Dealing with such 
matrices in the case of macroscopic surfaces [M ~ O(10 15 )] is not an easy task. But as 
we saw on the example of the creation and annihilation operators (jlj), the structure of 
the 2 M x 2 M matrix can be rather trivial consisting mostly from the Kronecker product 
of the 2x2 identity matrices. The Fock space formalism emphasizes the nontrivial 
aspects of the transition matrix leaving trivialities behind the scene. For example, a 
simple model of the surface growth consisting only in atomic deposition can be described 
within this formalism by the operator 

W F = FJ24, (12) 
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where F is the deposition rate. Using explicit representation for the state vectors (|2D 
and the operator identities ([7|) it is easy to see that the matrix element 

Wj = (a\W F |/3) = I ^ (13) 

is different from zero only provided |a) = aflP) for some site index % and the site 
is empty in state \/3) because according to (JTj) two atoms cannot sit on the same site. 
Thus, the transition matrix in this case is very sparse despite being very large. Similarly, 
hopping diffusion of a free atom on the surface can be described with the use of the 
operator 

W? = £>£>+ a,, (14) 

where D is the diffusion rate and we assumed that the atom that left site % may hop 
only to one of four nearest sites denoted here and in the following by the subscript N. 
Further, using the complete set of the system states, it is convenient to cast the 
ME (03J into the form 

d\t) 



where 



and 



dt T ^ < 15 > 

|t) = £P„(t)l<*) (16) 

a 



Wap-Safi^W. 



7a 
7 



(17) 



In practice one is usually interested in the average values of various densities described 
by the operators that are diagonal in the Fock space and can be calculated as 

Oit) = (6}(t) = J2 p a (t)(^\d\a). (18) 

a 

To express such averages in the form of matrix elements, as is conventional in quantum 
formalism, we introduce the matrix containing all possible states of the system with 
equal (unit) weight as 

|n) = J>>- (19) 

Now using the or t honor mality of the state vectors one gets 

0{t) = (Q\6\t). (20) 

Another important observation concerning vector |fi) in (TT9I) is that, as can be seen from 
( ITT]) , the sum of the matrix elements of matrix T over the left index is equal to zero: 
Yla Tap = 0. This is the consequence of the fundamental requirement of the probability 
conservation and in the quantum notation can be written as 

(fi|T = 0. (21) 
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Now using this identity and formally solving (j!5p via the operator exponent and 
assuming that the initial state at t — is the vacuum state (Q (the empty lattice) 
one can re-write (1201) as 



0(t) = (n\e- tT Oe tT \0). (22) 

From this we finally obtain an exact equation for the time evolution of the average value 
of 0(t): 

^ = (Q\[6,T]\t). (23) 

3. Evolution of island densities during submonolayer growth 

To derive the evolution equations for the island densities one needs to know the operators 
T and O entering (1231) . We first note that it suffices to know only the first non-diagonal 
term of T in (|T7|) because the density operators we are interested in are diagonal and 
thus will commute with the diagonal part in (T2"3"j) . (In case of necessity the diagonal part 
of T can be easily recovered from the non-diagonal part given below |33j.) The operator 
W that describes the epitaxial growth in the model under consideration consists of only 
two terms: the deposition term (112J1 and the diffusion term of the type of (1141) modified 
to include interatomic interactions through the dependence of the diffusion rate on the 
atom environment: 

This is done with the use of the activated dynamics that is conventionally used in the 
surface growth studies [MJ [15l [11] : 

A = Dexp[-(E N /kT) ^n % ] = D ]Jexp[-(E N /kT)n iN } 
= D J] (1 + [exp(-E N /kT) - 1] n lN ) En -^ T D]\n iN . (25) 

i N in 

Here En > is the energy of binding between neighbour atoms, in the second row 
use has been made of the property nf = rii, and n, = 1 — n,. The irreversible growth 
corresponds to the limit shown at the end of (|25l) so in the calculations below we will 
use the diffusion term in the form 

W D = DJ2<~^ (26) 

where according to (125J1 

a = cti]^[n iN . (27) 

Next we need to define the island density operators. We define atomic islands on 
the square lattice similar to the definition of the lattice animals [35], namely, as the 
clusters of atoms occupying the sites any pair of which can be connected by a path that 



Rate equations for epitaxial growth 



^e 



n 



^e 



i+y 



^e 



n.j 



Hi 



n 



^e 



i-\-x 



^e 



IT'i — y 



Figure 1. Encircled sites illustrate the structure of the monomer operator Ni t i in 
([^| . The crossed sites are the next neighbours (NN) to site i and the numbers are the 



values of rrii NN (the numbers of nucleation paths) in ([33 



traverses only nearest neighbour sites all belonging to the cluster. We characterize the 
island by its size s, configuration c, and the position on the lattice i. The latter can 
be defined arbitrarily but uniquely for all island of given configuration. The islands 
have different configurations if they cannot be obtained from each other by lattice 
translations. According to this definition, the dimers oriented along axis x and along 
axis y have different configurations. Thus, the density operator for the cluster of any size 
(including both the islands and the monomers) that has configuration c and is placed 
at cite i is 



N. 



(c) 



Uy n 



n 



jn- 



(28) 



{jec} {jNec} 
The second product in this expression is over the sites bordering the island by which 
we mean all sites that do not belong to the island but are nearest neighbour to at least 
one island site (we will denote the set of such sites as c). Operators rij N are needed 
to guarantee that the island contains exactly s atoms: without this factor the operator 

* (c) 

iV si could produce nonzero expectation value at configurations c that make parts of 
larger islands. In the case of monomer operator Nij the superscript c is not needed 
because the monomer configuration is unique (see figured]): 

N lti = m Yl n i+t> ( 29 ) 

e=±x,±y 

where x and y are the unit lattice vectors. 



3.1. Evolution equations for the island densities 

Substituting the island density operators (128]) into ( 123]) and calculating the commutator 
(see E]) one arrives on the evolution equations for their average values. The densities are 



Rate equations for epitaxial growth 9 

calculated according to the usual expressions: 

N^=M-^Nl% (30) 

i 

(N ltlNN N^) = M- 1 Y,(Ni,i N .N$), etc. (31) 

i 

In the first equation the carets over the symbols are omitted because here they are not 
operators but the number densities calculated according to (1201 ; in the second equation 
and in the following the subscripts i^ or «nn are considered to be relative coordinates 
with respect to the corresponding cluster or site (nearest- or next-nearest neighbours) 
and thus in ( 13T1) they are not averaged out. 

Thus, for the islands (s > 1) ( 123]) with the use of ( 1A.11J) and (1A.10I) takes the form 

^ = D J2' (N^N^}-Dj2mU^,m+Fj2 N { s ^-F7s c N^+. . . .(32) 

{c\l}i NN «NN {e\l} 

The prime over the first sum reminds that the summation over «nn is not over all NN 
sites to the s — 1-atom island but only over the NN positions that the detached atom can 
reach and remain free. By the dots we denoted the terms responsible for the coalescence. 
It would not be too difficult to write down these terms. Because of the symmetry of the 
square substrate a diffusing atom can be caught by one, two or three nearby clusters 
(islands and/or monomers). Similarly, the process of direct impingement can include 
up to four clusters. Thus, as can be seen, these processes could have been accounted 
for in ( T32l) with the use of 19 additional terms containing all admissible combinations of 
islands and monomers. Inclusion of these terms, however, would have made the equation 
very cumbersome. Besides, currently the research activity in the field of the irreversible 
growth is centred around the pre-coalescence regime where only the processes taken into 
account in ( l32|) participate. This regime corresponds to smaller coverages 9 < 0.2 and 
physically describes the ensemble of well separated isolated islands interspersed with 
low density monomer gas. Because the density of all clusters are very low the processes 
between closely spaced three or four clusters are small and can be neglected in this 
regime. If, however, there is need to account for such processes (at high coverages, for 
example), our technique allows for rigorous derivation of necessary equations. But in 
the present paper we will neglect the coalescence terms in all equations below. 

In order to include spatial correlations and/or fluctuations into the mean-field 
description, equation describing the monomer diffusion have been proposed [I5j [23j 
36] . To obtain similar equation in our approach we postpone the site averaging in 
301) corresponding to Nij till later. Now from ( 1 A. 1411 and ( 1A.12I) one gets 



FN 0>i - AFN ltl + D I £) V Mn - 4iV M )~J2 
\ in / 'nn 



^ = FN 0)i - 4FiV lj4 + D ( > ] N 1)is - 4iV M | - > ] m« NN (iV MNN iV M ), (33) 



where A*o,j is defined in ( 1A.13I) and the prime on the summation over zn means that the 
corresponding terms are present only if site i does not have nearest neighbour atoms. 
This is a consequence of the irreversible dynamics we are studying. If the atom at % had 
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a neighbour it would have been bound and so could not produce a free monomer via 
the diffusion hop. 

The comparison of ( |33l) with the equations proposed in [T5(, |28] as well as derivation 
with its help of the conventional site averaged rate equation for Ni will be given in the 
next section. 

4. Rate equations 

The equations (132]) and (133]) cannot be used directly for the calculation of island densities 
because they contain unknown terms — the averages of the operator products in the 
angular brackets. The standard way to proceed in such cases is to decouple them by 
replacing the mean of the operator product with the product of the means 

(Ni^N^) « JViJVjO. (34) 

This is the conventional mean field approximation that proved its efficiency in a lot of 
many-body problems. In the theory of epitaxial growth, however, the standard mean- 
field approach was found to be too crude [2]. Therefore, it was suggested to account for 
the correlations neglected in ( 134]) via the use of parameters o~s defined according to the 
equality 

T l {Ni tim Ni e) ) = <ri s) N 1 NJrl (35) 

Because of their physical meaning the parameters are called the capture numbers [2]. 
Substituting ( 135]) into ( 1321) and augmenting the set by the equation for the monomer 
density derived in 14.11 we arrive at a formally closed set of equations that describe the 
evolution of the island densities. The set, however, is not complete without explicit 
expressions for the capture numbers. As we noted in the Introduction, there exist many 
heuristic approaches aiming to define the CNs independently. Our aim in the present 
paper is to propose a method to derive the exact values of the CNs from the KMC 
simulation data. This would provide us with a definite mean of assessing the quality 
of the CNs obtained in empirical approaches. We achieve this by noting that equations 
( 132]) and (133]) are rigorously derived from the ME (ITTj) that exactly describes the KMC 
stochastic process. Thus, the KMC data should satisfy these equations exactly hence 
by comparing them with the empirical REs one can calculate the exact values of CNs. 
For example, (135]) can be trivially solved as 

^ = ^(NiAnNtoyKN®. (36) 

»NN 

We note that all quantities on the right hand side of this expression can be measured 
in the course of exact KMC simulations thus providing the exact values of CNs. A 
comment is in order concerning the influence on this result of the coalescence terms 
that we omitted in our evolution equations. It should be noted that they also are 
discarded in the REs that we are studying in this paper [5] [2] and that the CNs cannot 
be calculated with better accuracy than their definition allows. Thus, for the equations 
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derived in the "no coalescence" approximation in the KMC simulations we may count 
only those configurations where exactly one monomer and one island participate and 
discard those processes which lead to coalescence. In this sense the measurements of 
CNs with the use of (1361) and KMC may be considered to be in exact correspondence 
with its definition. 

Individual capture numbers for every island configuration introduced in fl35|) and 
( T36|) can be efficient in the cases when the configurations are easily enumerable, e. g., in 
the case of one- dimensional islands either on the one-dimensional (ID) substrate or on 
2D substrates in the cases of very anisotropic quasi-lD growth [371 E3 H]- Because ID 
islands have only one configuration characterized by their size s, the superscript (c) for 
such islands in the above equations is superfluous. 

In more common case of isotropic growth in 2D, however, the number of island 
configurations is so large that for s > 50 even their enumeration is infeasible [35]. 
In such cases cruder description of island morphologies is needed. Usually only the 
island size statistic is being gathered both in the KMC simulations and in the growth 
experiments. So now we derive the formulas generalizing (136]) to this case. 

4-1. REs in isotropic case 

Because the number of configurations of islands of sizes s > 50 is unmanageable [35], 
the rate equations are usually written for the total density of islands containing the 
same number of atoms as 

iV s = ^iVt). (37) 

c 

The standard REs for the island densities during the irreversible growth have the form 

m 

— -i = DN x (<r,-iN,-i - a s N s ) + Fk.-xN...i - Fk s N S} s = 2, 3, . . . (38) 

We can derive them from (132]) by summing the latter over the configurations 
corresponding to the same island size as in ([37]) . In this way we obtain the left hand 
side of (J38]) . But on the right hand side one finds the terms of the kind 

-Fj>iVt) (39) 

c 

[see the last term in (132]) ] Identifying it with the last term in (138]) we arrive at the 
following definition of the capture numbers for direct impingement 

K s = J2$cNi c) /N s . (40) 

c 

One may wonder whether this definition is consistent with the positive deposition terms 
in (]32]) and f]38]) because the positive term in ( ]32]) has the form quite different from that 
of the negative term. To see that they acquire the same form after the summation over 
configurations it is sufficient to note that every term of the kind AQl\ will enter the 
sum proportional to F in (1321) as many times as there is configurations with s atoms 
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that would lead to the configuration (c \ 1) when the atom is taken away. Obviously 
that the atoms can be taken only from one of the sites in the configuration border so 
there is exactly (s — l)( c \i) such terms. 

Similarly, using the second term on the first line of ( 132]) and comparing it with (|38j) 
one gets 

«* = E m ^ (fti^fiU/KN.. (41) 

«NN,c 

The contribution of the first (positive) term in ( 132]) can be analysed similar to the direct 
impingement case discussed above. 

To close the set of the rate equations for the island densities, the equation for 
Ni = M _1 J2i -Wi,t is needed. Its left hand side is trivially obtained by averaging the left 
hand side of (133]) over the surface, so let us consider in detail the right hand side of this 
equation. The first term the "empty island" operator gives a nonzero contribution equal 
to unity only when site i and all its neighbours are empty. This excludes all occupied 
sites as well as the sites bordering the islands and the monomers. Taking into account 
(j4"0|) and the second term in (1331) which is simply an additional monomer contribution 
— k±Ni («i = 4) we obtain the term in the parentheses of the conventional RE for the 
monomers (see, e. g., [25]): 

<^± = F ( i _ _ 2/Cl JVi - > > s iV s ) - 2DatN'f - DN,} ] a s N s . (42) 



2/tiiVi - Y^ K * N s ) ~ ZDaxNl - DNi J2 a s N s- 

S>1 J S>1 



The diffusion term in ( 133]) is the discrete Laplacian that consists of two terms: the 
positive contribution from the monomer hops on site i from four neighbour sites with 
monomer densities iVi ( i N and the negative contribution from site % on the neighbour sites. 
Thus, when summed over the whole lattice the negative terms tend to be cancelled by 
the positive contributions from the neighbour sites. In the presence of islands and other 
monomers, however, the positive terms may not exist because of the occupied sites 
nearest to them. This leaves uncompensated terms in the negative contribution that 
equals exactly to the number m iNN of the nucleation paths to the mentioned islands or 
monomers [23]. This in accordance with (141]) gives the last two negative terms in ( 142]) . 
Finally, the contribution that gives factor 2 to <j\ term comes from the last term of (13"3"]) . 

4-2. Comparison with equations from the self- consistent approach to the CNs 
calculation 115] 



To facilitate the comparison we first rewrite ( 133]) in the continuous substrate coordinates 
i — > x as 

dN ^ ~ FiVo(x) - KiFNxix) + DV 2 iVi(x) - ZMiV^x) 2 , (43) 

where iVo(x) is the continuous version of the "empty island" operator Nqj from ( 133] ). 
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In this notation equation (11) from [15J would read (we rearranged some terms for 
convenience) 



dNx( 



dt 



F - KiFNxix) + DV^^x) - 2Da 1 N x N 1 (x) 

oo oo 

- FJ2^N S - DNiW^tTsNs, (44) 



S>1 



where 



^X^/^M^IX' (45) 



(the lattice constant a is assumed to be equal to unity). In comparing (144"]) with (143]) 
we first note that the last term in (1431) that describes the island nucleation is local in x. 
This is reasonable because the adatoms nucleate an island when meeting at one spatial 
point in the continuum approximation (or at the nearest sites on the lattice) while 
according to ( 1441) and f j4"5]) the adatom at point x nucleates an island with any other 
adatom in the system which is unphysical. Secondly, in ( 1441) the nucleation term enters 
with the coefficient 2, in contrast to our equation ( 143]) . The origin of this discrepancy is 
connected with the last term of (1441) that is absent in our equation (1431 but is present as 
the last term in (142]) ). As we saw in section I4TT1 it arises from the averaging the discrete 
Laplacian in ( 1331) over the whole surface. Thus, in the continuum approximation this 
term also should appear from the integration of V 2 iVi(x) = div [V-/V\(x)] as the sink 
term for the monomer attachment to the island boundaries. But as it is already present 
in (1441) this means that the capture of monomers is counted twice in this equation, unless 
some very special treatment of the boundary conditions is undertaken which seems to 
be the case in [15]. 

Similarly, the first term on the second line in (1441) appears in our formalism only 
after the first term in ( 1431) is averaged over the whole surface, as in (142]) ). These direct 
impingement terms are less important under the usual growth conditions and are often 
neglected, though they may became important in some cases. Apart from these terms, 
our equation (1431) is very similar to the % = 1 case of the equation derived in [271128] . In 
this approach the nucleation is local and the capture term DNi(x) Y^>i °~sN s is absent 
from equation (2.7) in [28]. The nucleation term, however, enters with the coefficient 
i + 1 = 2 which we consider to be double counting. Careful comparison of all three 
equations with exact KMC simulations is necessary to clarify this complicated issue. 

5. The nucleation rate 

To illustrate the above formalism, in this section we present the calculation of <J\ with 
the use of the KMC simulations and the formula 

*i= E "^.(VmW) ( 46 ) 

jnngnn 



Rate equations for epitaxial growth 14 

that can be derived from the last term in (1331) . as explained in 14.11 As was noted 
in the Introduction, <7\ is the most important of the CNs in the systems belonging 
to the so-called i=l case when all islands with s > 1 are stable. This is because in 
the usually considered case of large R when the direct impingement terms are small, 
<7i effectively defines the nucleation of dimers and hence the total number of islands 
in the system because every island nucleates as a dimer in the i = 1 case [15] . The 
total island density measured experimentally can be used to define the diffusivity of 
atomic monomers on the surface [ID] . Therefore, an accurate theoretical calculation of 
o~\ could lead to more reliable empirical estimates of the microscopic parameters that 
define the surface diffusivity. Because of this, <Ti was calculated by several authors but 
controversies remained. For example, the widely cited analytic expression [391 13 HSl 12] 

47T 

ai * hcrnTy (47) 

where C is a numerical constant. In the steady-state regime at large R the monomer 
density behaves as [5] 

Nt ~ (3n 2 6R 2 )- l/3 . (48) 

It should be noted that this formula is non-physical at small coverage where N\ ~ 6 
while according to (148]) N\ diverges. Therefore, we will use it only above the crossover 
coverage separating the transient and the steady-state regimes [2] 

6* = R~ 1/2 (49) 

(see figure [2]). Thus, substituting (J4"8"]l into (J4TJ) one gets 

ai ~ \nC 1 + (\nR-\n9)/3 , (50) 

where C\ = C/(37r 2 ) 1 / 3 . As is seen, the dependence on R in (150]) is formally as strong as 
on 9 and in practice is even stronger because in experiments and simulations on growth 
in precoalescence regime the coverage usually varies within one order of magnitude 
(0 ~ 1 — 20%) while R varies over several orders of magnitude. For example, in 
simulations of [25J) R = 10 5 — 10 8 . Yet the authors arrive at the conclusion that 
the coverage dependence is more important. 

To compare the values of <J\ calculated in the framework of the above theories with 
the exact values given by (|46|) we performed KMC simulations at two values of R and 
at small coverage shown in figures [2] and [3j The choice of parameters was dictated by 
the necessity to obtain reasonable statistics which depend on the monomer density Ni, 
as can be seen from ( |46J) . The latter diminishes with the growth of both R and 6, as is 
seen from ( 14"8"J) . Therefore, we have chosen the region near the crossover coverage (T4"9"|) 
where the density is at its maximum (see figure |2]) [2]. 

In the harder R = 10 5 case, 400 KMC runs in the system with 4096 x 4096 sites 
were made. In the easier R = 10 4 case 2048 x 2048 systems were used and lesser 
statistics were gathered to obtain similar statistical errors. No finite size effects were 
noticed. As can be seen from figure [3j ( T4T1) with the value C = 1 taken from [15] 
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Figure 2. Thick lines: the monomer densities obtained in KMC simulations for two 
values of the diffusion to deposition rates ratio R = 10 4 (upper line) and R — 10 5 (lower 
line); thin solid lines: the densities in the steady-state regime calculated according to 
the asymptotic formula (|50|) : vertical dashed lines show the positions of the crossover 
coverage 9* from (j49|l . 
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Figure 3. Filled symbols: KMC measurements with the use of the exact formula (|23|) 
at diffusion to deposition rates ratios R — 10 4 (•) and R — 10 5 (■); similar empty 
symbols: KMC measurements with the use of the procedure described in [5D] and [25] , 
Thin solid lines: <j\ calculated according to (H7|) with C = 1 [TS] with the use of the 
simulated values of Ni (see figure [2]); dashed lines: ([501) for 9 > 9* . Thick solid line: 
fit to the R = 10 7 KMC data from figure 2 of (25]. 
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semiquantitatively describes the simulated data over the whole range of coverage studied 
both in the transient regime and in the course of the steady-state growth, especially at 
the larger diffusion to deposition rates ratio R, as predicted theoretically. In the paper 
[39] similar value of C ~ 1.38 was found. Thus, the results of [391 US] an d our results 
seems to rule out the values C ~ 32 — 64 advocated in [2E1 E]. 

On the other hand, the KMC procedure suggested in [20] and [25] produces results in 
qualitative disagreement with the exact data in the transient region, though practically 
coincides with them within the statistical errors in the steady-state regime, as can be seen 
from figure [3j The reason for this seems to be clear. The procedure in question consists 
in stopping the deposition from the external source at any given moment and detaching 
the monomers captured by the islands or separating the nucleated dimer and putting 
them back into the homogeneous deposition flux over the free sites without occupied 
neighbours [201 ES] . The capture/nucleation events are counted for every island size. By 
continuing this arrested growth simulation sufficiently long one can obtain arbitrarily 
good statistics on capture numbers. Obviously, however, that this procedure may give 
quite accurate results only in the steady state regime where the diffusion profile iVi(x) 
is time-independent, so the origin of the homogeneously deposited atoms, — either from 
the physical deposition flux or from the artificially simulated one, — is irrelevant: the 
profile is always the same |17j . 

At small coverages at the early stage of the growth, however, the diffusion profile 
only strives to reach the steady-state shape and so is time dependent. Thus, at the start 
of the arrested deposition the diffusion profile has a physical shape which permanently 
changes with time and as the arrested deposition simulation continues the shape of the 
profile in the absence of the external deposition became unphysical and the simulation 
produces incorrect G\ values, as can be seen from figure [3] 

6. Conclusion 

In this paper we used the second quantization representation of the stochastic master 
equation governing the coherent epitaxial growth to develop a formalism that allows for 
a rigorous derivation of the rate equations for the epitaxial growth. The REs obtained in 
sec.[3]are more detailed than those found in literature in that they describe the evolution 
of individual island configurations. This may be of value in the cases when the number 
of configurations is restricted, as in the case of ID growth [3?! [381 II]- I n the majority 
of cases, however, the number of the island configurations is too large to be used in 
practical calculations [35]. Therefore, it is conventional to use the REs where the island 
densities comprising all island configurations of the same size are used. We illustrated 
the derivation of such REs in section 14.11 In some cases, however, experimental data 
show that not all configurations are equal. The so-called magic clusters may play more 
significant role than islands of other shapes [4"0~ 1 HI ] . Our equations of Sec. [3] are capable 
of assuring a special treatment to such configurations. Furthermore, in the case of 
quantum dot growth, islands of different morphologies exhibit qualitatively different 
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growth behaviours [421 S3] which also would require more detailed characterization of 
islands than only by their size. Finally, the detailed derivation of the REs can be 
helpful in better understanding the approximations underlying the REs and the means 
of improving them with the use of more accurate decouplings. 

As we showed in the present paper, even in the simplest case of the irreversible 
submonolayer growth in precoalescence regime that has been extensively investigated for 
about forty years there still remain controversial issues that need further investigation. 
Obviously that in more complex cases such as the currently hot subject of the quantum 
dot growth [121 H3] the difficulties will enlarge manifold. We expect that the formalism 
developed in the present paper will prove to be an indispensable tool of derivation of 
the REs for such cases. The derivation in the framework of the second quantisation 
formalism, though cumbersome, is quite straightforward. Besides, it is very flexible 
and do not reduces to the simple case of the hard core bosons we considered in the 
present paper. Multiple site occupation, coalescence, 3D growth, reversibility, etc all 
can be accounted for both in the quantum formalism [291 ES EH [321 [33] and in the REs 



Appendix 

The commutators in (1231) can be calculated with the use of the chain rule 

[A X A 2 ...A k ,X]=[A 1 ,X]A 2 ...A k + A x [A 2 ,X]A 3 ...A k + ... 

+ A 1 ...A k _ 1 [A k ,X], (A.l) 

the definitions 

rii = afai, fii = 1 — rii = a^af (A. 2) 

and the commutators 

[m, of] = - [fii, at] = at (A.3) 

[rii, a t ] = - [fii, Oi] = -a*. (A.4) 

As can be easily seen, in the case of binary site occupancy vector \fi) in (TT9l can 
be expressed through vectors ([3]) as 

|n> = n(l°>< + l 1 >0- (A.5) 



From this representation it is easy to see that inside the matrix elements of type ( 12011 . 
i. e., with vector \Q) on the left, the nondiagonal operators can always be turned into 
diagonal ones as 

(Q\at = (n\fii (A.6) 

(n\oi = (n\m. (A.7) 

Care should be taken to keep noncommuting operators in correct order. For example, 
vital to the results below is that operator a* in ( TT7I) is the rightmost in the product. The 
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evolution equations that we derive strongly depend on this because of the following not 
quite intuitive identities: 

(fl\fiiai = (Q\rii (A. 8) 

(ttlfnaf = (A.9) 

[cf. flAJl— dA21)]. 



Because physically there exists a qualitative difference between the islands and 
mobile monomers, it should be reflected in the formalism. Therefore, we will discuss 
the calculation of the commutators for these entities separately. 

Appendix A.l. The island commutators (s > 1) 

First we calculate the simpler commutator of the deposition term in ( |2^|) with the island 
operator (1281) . From ( 1A.3I) it is seen that the commutator of the creation operators with 
the island operator (128]) will consist of two parts : the commutator with operators n^ on 
interior island sites and with operators ni on the border sites. Though the commutators 
in ( 1A.3I) differ only by sign, the results are somewhat different. The border commutators 



when acting on the vector (Q\ result in the same expression as the island operator ( j28|) 
itself while the interior commutators lead to islands with one less atom that arise when 
every n^ is sequentially replaced with fii corresponding to an empty site. We denote 
symbolically the set of such islands as {c \ 1}: 

(n\[N$ , W F ] = F(Q\ I J2 AtYj - sN^ ) , (A.10) 

\{c\l} / 

where the first sum is over s islands whose configurations differ from c by the absence 
of each of s atoms constituting the initial island. It is to be noted that in some cases 
this might lead to disconnected sets, i. e., to the appearance of two or more separate 
islands or monomers. Such cases correspond to the growth with the coalescence and we 
do not consider them in the present paper devoted to pre-coalescence growth regime. 
Their inclusion would lead to exact equations of more general type. 

The calculation of the commutator with the second term in ( 12^1) in the case of 
islands is very similar due to the structure of fij in fl26|) . Because the annihilation 
operator in W D is multiplied by factors fii NN on all neighbour sites, fij will commute 
with any diagonal operator on site i that contains at least one of n« NN as a factor. But 
such are all operators entering N^ c - by definition. Therefore, in the commutator cij can 
be treated as a number similar to F and the commutator calculated exactly as in (1A.10I) . 
The main complication is connected with acting with D^V a iN on the vector (Q\, But 
as explained above, for every boundary site only the sites that have no neighbours with 
Hi will contribute, that is, only the NN sites to the island with configuration c and the 
NN sites to the islands with one less atom [see (1A.10J) ] amounting to 



{c\l}i NN »NN 



Rate equations for epitaxial growth 19 

where the prime over the first sum is to indicate that the summation over the NN sites 
comprises only those sites (there may be from one to three) which the detached atom 
can rich in one hop and rrii NN is the number of nucleation paths for the NN atom to 
reach the island 



Appendix A. 2. The monomer case 

Because the monomer is a limiting case of an island, we can proceed with the calculation 
of commutator as above by only bearing in mind that monomers differ from islands in 
two respects. First, when an atom is taken away from the monomer in the first terms 
in ( 1A.10J) and ( 1A.11J) we obtain an "empty island" . In the first case we find 



(n\[N 1>h W F ] = F(n\ (N ,i - 4A> M ) , (A.12) 

where the "empty island" operator 

N ,i = fii Y\ n i+i (A. 13) 

e=±x,±y 

has formal structure of N\ t i ( 1291) shown in figure Q] but with the central operator n^ 
replaced with n^. The second commutator takes the form 

(n\[N hh W D ] = (n\(-WN lti ) + (Q\D (j^^^f - ^m JNN iV MNN iV M j , (A.14) 

\ IN 'NN / 

where the last two terms in parentheses are analogous to the terms entering ( 1A.11I) for 



the islands while the first term is due to the commutator with a from W D (}26l) that is 
absent in the island case. The operator Nq f is similar to the empty island operator Nq^ 
and can be obtained from it by means of (IA.8J) . namely by replacing one of fii N factors 
surrounding site i with n iN . Its action is to suppress those configurations in which an 
atom detaches from an island. 
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